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ABSTRACT 

We present a near-infrared extinction study of nine dense cores at evolutionary 
stages between starless to Class I. Our results show that the density structure 
of all but one observed cores can be modeled with a single power law p oc r p 
between ~ 0.2R — R of the cores. The starless cores in our sample show two 
different types of density structures, one follows p ~ —1.0 and the other follows 
p ~ —2.5, while the protostellar cores all have p ~ —2.5. The similarity between 
the prestellar cores with p ~ —2.5 and protostellar cores implies that those 
prestellar cores could be evolving towards the protostellar stage. The slope of 
p ~ —2.5 is steeper than that of an singular isothermal sphere, which may be 
interpreted with the evolutionary model of cores with finite mass. 

Subject headings: dust, extinction — ISM: globules — ISM: clouds — stars: formation 
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1. Introduction 

To understand the nature of star formation processes, it is crucial to measure the 
fundamental physical parameters of molecular cores during their evolution. Density 
structure is one of the most important parameters which can provide critical information 
to test theoretical models. The simplest model for a pre-collapsing core is a Bonnor-Ebert 
sphere (hereafter BES; Bonnor 1956 and Ebert 1955), which describes an isothermal, 
non-magnetized, non-rotating, hydrostatic core bounded by external pressure. For a 
collapsing core, the standard model is inside-out collapse model which predicts the density 
p oc r~ 2 before the gravitational collapse sets in and p oc r~ L5 in the inner free fall region 
(Shu 1977). Vorobyov and Basu (2005) find the power-law index of density profile becomes 
much steeper than p oc r~ 2 in the outer region of the core if the core has finite mass 
reservoir. 

Two methods have been commonly used to map the density structure of cores. One 
is to measure the thermal dust emission at millimeter or submillimeter wavelengths using 
single dishes or interferometers, and the other is to measure the extinction of the background 
stars at near infrared (NIR) wavelengths. Large surveys of core density structure can be 
conducted relatively easily using single dishes, but the resolutions are often limited and the 
density profile can depend on observing wavelength (Shirley et al. 2000; Kauffmann et al. 
2008). Interferometers can be used to achieve high resolution in the inner-most region, but 
the overall structure are often resolved out (Harvey et al. 2003a). Given enough background 
stars and integration time, NIR extinction observations have the best potential to achieve 
both high resolution and large scale mapping of the core density structure (Alves et al. 
2001). 

Previous studies have shown that BES or power-law density profiles are good 
descriptions for density profiles of molecular cores. Kandori et al. (2005) found that Bok 
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Globlules can be described as either gravitationally stable or unstable BESs with NIR 
extinction observations. Shirley et al. (2000) found that the density profiles of 21 low-mass 
cores from starless to Class I stages can be fitted with power laws p oc r~ 15 to p oc r~ 2S in 
the resolved regions with assumed dust temperature Td{r) oc r~ 0A . It is unclear, however, 
whether or not the variations in the observed density structure are due to evolutionary 
effects. In this study, we focus on investigating the evolutionary trend of core density 
profiles with NIR extinction observations. The observations and data reduction processes 
are described in §2. The derived visual extinction {Ay) distribution and the radial density 
profiles are presented in §3. We discuss the implications from our results in §4. 

2. Observation and Data Reduction 

We observed five isolated cores (CB68, L483, CB188, L158 and B59) and two fields 
in Serpens Molecular Cloud in the H and Ks bands using the Wide-field InfraRed Camera 
(WIRCam) on Canada-France-Hawaii Telescope (CFHT) in 2008 April and May. These 
fields are located in low Galactic latitude and hence they have abundant background stars. 
The isolated cores and the selected dense cores in two Serpens fields have relatively round 
shape and contains only a protostar (Class or I), except for B59. B59 is a cluster forming 
core and is chosen for the comparison between single star and cluster formation. 

For each observing field, we performed in total a 17.5 minute integration (5 sets 
of 21 dithered images) in H band and a 3.5 minute integration (one set of 21 dithered 
images) in Ks band. Standard bias subtraction, flat-fielding, bad pixel masking, and sky 
background subtraction were performed using Tiwi system, the WIRCam pipeline Q The 
sky background of each image was subtracted based on the co-moving averaged frame with 
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length of 15 minutes. We combined all images of each observing field and used SExtractor 
package (Bertin and Arnouts 1996) for the source extraction and photometry except for the 
extremly crowded CB188 field for which we used DAOPHOT package (Stetson 1987). The 
photometry results from both methods were calibrated using the point sources with 2MASS 
detections between 13 to 15 magnitude. The typical limiting magnitudes with 3 a detection 
are about 20.7 and 19.5 in H and Ks band. 

3. Analysis and Results 
3.1. The Visual Extinction Ay 

The Ay value is determined from near-infrared color excess Eh-ks of stars in the H 
and Ks bands (NICE method, Lada et al. 1994), i.e., Ay = r^' Ks E H _ Ks . r ^~ Ks is given by 
a dust reddening law and Eh-ks is evaluated from (H — Ks) b se rved — (H — Ks)i ntr insic- We 
use r^~ Ks = 19.245 as derived from the Rv=5.5 model in Weingartner and Draine (2001, 
WD). Recent studies suggest that the reddening law in dense cores is better described 
with WD's Rv=5.5 model rather than their Rv=3.1 model which fits the reddening law 
derived from more diffused regions (Indebetouw et al. 2005; Chapman et al. 2009). E H -k s 
can be obtained from observations once (H — Ks)i ntr insic of the star is known. We adopt 
(H — Ks) intr insic = 0.13 mag (Lada et al. 1994) and assume it is a constant for all observed 
fields. 

The detected point sources in our fields are not necessary stars; they could be young 
stellar objects or distant galaxies. We cannot perform the standard JHK S color-color 
diagram to select stars since we did not observe our target sources in J band to save 
observation time. We can however take advantage of the existing catalogs from the Spitzer 
c2d Legacy Project (Evans et al. 2003) to select bona-fide stars detected in our fields. 
The c2d Project has developed a sophisticated source classification procedure based on the 
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2MASS and Spitzer data (Harvey et al. 2007). If the point sources are not contained in c2d 
catalog, we assume them to be stars. The derived extinction maps are shown in Figure 1. 

3.2. The Density Profile 

We assume the three dimensional shape of our observing cores are spheres and derive 
their density profiles. We first determine the center, radius, and background extinction of 
the cores (Table 1), then annularly average the derived extinction to produce the extinction 
profile (Figure 2). 

The centers of the four cores with each a single embedded star (CB68, L483, CB188 
and Serp-Bolo20) are set to be the position of their associated protostar since the protostar 
and the peak of Ay distribution almost coincide. The center of the four starless cores 
(L158, Serp-Bolol, Serp-Bolo5 and Serp-Bolol9) and the cluster forming core, B59, is set at 
the peak of their Ay distribution. We set the annular averaging step as Ar and ensure that 
each step contains at least five detections to obtain properly the average Ay. Foreground 
star contamination may affect the derived extinction. We expect, however, the foreground 
star fraction is small since all the selected regions are nearby (~ 200 pc). If foreground 
stars contaminate the data in Ar, the derived Ay of foreground stars will deviate from 
the median of the data. Therefore, to remove the possible bias from foreground stars, 
we compute the median (m^ v ) and the standard deviation (oa v ) of Ay distribution in 
each annulus and remove the sources with Ay larger than rriA v + 1.5 <ja v or smaller than 
rriA v — 1.5 cr Av . Then we annularly average Ay in each Ar. 

At a certain distance from the center, the Ay approaches nearly a constant value which 
is defined as the background Ay. We calculate the average background Ay (BG) and the 
standard deviation (SD) in the background region. The radius of the cores is defined as the 
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distance from center to the position where the Ay exceeds BG + 3SD and the Ay profiles 
are obtained by subtracting BG. The relationship between the Ay and H 2 column density 
are derived from (JV(ifj) + 2N(H 2 ))/E B -v = 5.8 x 10 21 cm" 2 mag" 1 (Bohlin et al. 1978). 
Assuming all hydrogen in the dense cores is in H 2 form and adopting the Ry=5.5 reddening 
law (where R v = A V /E B _ V ), we obtain N(H 2 )/A V = 5.27 x 10 20 cm" 2 mag" 1 . 



The observed column density profiles can be modeled with various functions that 
represent different inherent physical conditions. Bonnor-Ebert spheres are commonly used 
for presumed hydrostatic cores with flat density profiles in the center and are truncated at 
radius R. On the other hand, the solutions from different collapsing models could be as 
simple as a single power-law function or as complex as non-analytic functions. In general, 
these solutions can have approximate forms described by a power-law function with different 
power-law indices at different positions from the center. In this work, we model the density 
structure of dense cores with a single power-law profile as the first approach and further 
examine the power-law index variation within the density profile in §4.2. 

Assuming the volume density p{r) oc r p , the observed column density N(w) at the 
projected radius w in the plane of the sky is 



The best fit power-law index p in p{r) is obtained by minimizing the \ 2 between derived 
density profile and power-law density model. We adopt the Bootstrap test to determine the 
uncertainty in p (see astronomical application in Wallin et al. 2007 and "Numerical Recipes" 
by Press et al. 1993). The main idea of the Bootstrap test is to find the uncertainty of any 
operation as the dispersion of the results from the same operation on a large number of 
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simulated data sets. We produce 100 sets of simulated density profiles for each core, where 
each data point in the simulated profiles is generated from a Gaussian random number 
generator with the mean and standard deviation of the probability distribution equal to 
the mean and the uncertainty of the observed data. We therefore obtain the uncertainty 
in p from performing our fitting procedure on each simulated data set and measuring the 
standard deviation of the fitting results for each core. The best-fit results are listed in Table 
1. Figure 2 shows the best-fit results overploted on the observed data. 

Most of our samples have data in the range of 0.2R — R except for Serp-Bolol 
(~ 0.3R — R) and CB188 (~ 0.06i? — R). To compare fairly the average power-law index 
of each core with the same fitting range, we further fit all samples in 0.3R — R and fit 
CB188 in 0.2R — R and 0.3R — R. Figure 3 shows the original best-fit results and the 
results obtained with the fitting range in 0.2R — R and 0.3R — R. The normalized x square, 
X%, of the fitting results are all smaller than one because the dispersion of the data in the 
same radius are in general larger than the fitting uncertainty. It is common to use the 
measurement uncertainty as the weighting in fitting processes in many works (e.g. Harvey 
et al. 2003b). We use however the dispersion of annularly averaged data as the error bar 
here to include the uncertainty in assuming spherical structure. This assumption reduces 
the fitting x%, but increases the uncertainty in p. 

4. Discussion 
4.1. The evolution of the density structure 

We investigate the power-law index variation of the density structure at different 
evolutionary stages. To increase the size of our sample, we include the results from several 
similar infrared extinction studies, including L694-2 (Harvey et al. 2003b), B335 (Harvey 
et al. 2001), Thumbprint Nebula and DC303.8 (Kainulainen et al. 2007). In these works, 
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the extinction values are all determined from H-K color excess and the density profiles are 
all fitted with single power law, which is the same as our approach. Although the computer 
codes used to do the fitting may not be identical, we expect that these differences will not 
be significant. 

The six starless cores in this sample show two different types of density structure. 
Serp-Bolol9 and L694-2 have the best-fit power-law index ~ -2.5 (hereafter "steep type"), 
while other four cores (Serp-Bolo5, Serp-Bolol, L158 and Thumbprint Nebula) have the 
power-law index ~ -f.O to -f.4 (hereafter "shallow type"). This index range in shallow 
type reduces to ~ -f.O to -f .2 if we reduce the fitting range to 0.2R — R. The difference 
between these two types is significant even after considering the uncertainty of best-fit 
results and it is obtained under fair comparison with similar fitting range. On the other 
hand, the power-law indices of Class samples do not have clear two different types as 
starless cores but a range about -f.9 to -2.6. The two samples in Class I stage show two 
different power-law indices again (-2.7 in Serp-Bolo20 and -f .2 in CBf88). This discrepancy 
may result from different fitting ranges. Indeed, their best-fit results become similar to 
Class sources after considering only a fixed fitting range within ~ 0.2R — R (Figure 3). 

For the starless cores, the Ay difference between the center and the edge also show 
a dichotomous nature. The shallow type starless cores have much smaller Ay difference 
(~ 22 mag) compared to the steep type starless and protostellar cores (~ 45 mag). If all 
starless cores in our sample eventually evolve into protostellar cores, our results suggest 
that the starless cores with shallower density structure will evolve into steeper ones via 
certain core-forming processes and their density structure have a averaged p ~ —2.5 during 
the formation of the protostars. One way to examine if the starless cores will evolve into 
the protostellar cores is to probe their kinematic structure. Lee et al. (2004) shows that 
both L158 (shallow type) and L694-2 (steep type) have infall signatures in the molecular 
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spectra, which are consistent with our suggestions that some shallow type starless cores will 
evolve to the steep type during the formation of protostars. 

4.2. Power-law index variation within density profile 

Some core collapsing models predict a power-law index variation within the density 
profile. For example, a discontinuity in power-law indices of the density profile separates the 
interior infalling region from the exterior static region, and the discontinuity will propagate 
outward with sound speed during collapse (e.g. Shu 1977, Fatuzzo et al. 2004). Since the 
crossing time of the discontinuity (~ 10 5 yr, estimated with T = 10K, m = 2.33 amu and 
R = 0.1 pc) is comparable to the statistical lifetime of Class and Class I stages 10 4 yr 
and 10 5 yr, Andre et al. 2000), we should be able to see the discontinuity in the observed 
density profiles in Class or Class I stages. 

The best-fit results with fitting range of 0.3R — R showed in Figure 3 are in general 
steeper than fitting range of 0.2R — R, suggesting the power-law index may vary within 
core density profile and specifically, it may become steeper near the edge. We examine 
this trend by modified fitting processes. Considering the simplest case, a core has a two 
power-law volume density profile with the discontinuity at r = R'. Its column density 
profile within projected radius R' will be an integration of the volume density over regions 
with both power-law indices along the line of sight while the column density outside of 
R' will only reflect the power-law index in the outer region. Therefore, we can obtain the 
best-fit power-law index p 1 in the region Rl < r < R and then use pi as a known parameter 
to determine the best-fit power-law index p 2 in the r < R' region. R' is the radius chosen 
where the single power-law fitting in the region of R! — R reaches minimum Xn- 

The feasibility of the above procedures, however, is limited due to the large uncertainty 
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of annularly averaged data. When we vary R' from outer region toward inner region, we 
may not be able to find obvious minimum Xn an d the best fit p\ may vary significantly. If 
we set the criteria for variation of best-fit p± smaller than ±0.3 (i.e., the typical uncertainty 
in best- fit single power-law profile), L483 is the only possible candidate that shows deviation 
from a single power-law profile. In L483, the Xn °f the two power-law model are ~ six times 
smaller than the x% °f single power-law profile. This difference, however, is not significant 
since the best fit results of two models are both obtained with x"n < 1- The best-fit two 
power-law profile of L483 has the steep power index pi ~ —3.9 at r > Rl (~ 0.5R) and 
p 2 ~ — 1.7 at r < R. This result shows the possible power-law index variation in a Class 
sample, the power-law index —1.7 in inner region is consistent with the power-law index in 
the infall region of inside-out collapsing scenario. The power-law index, however, becomes 
steeper than — 2 at r > 0.5R suggesting other effects (e.g. see the discussions in §4.3) 
besides inside-out collapse should be considered in outer parts of the core. 

4.3. The implication for core evolution 

Our results show that the density profiles evolve from a shallower one with p ~ — 1.0 
to a steeper one with p ~ —2.5 in the epoch of core formation and the steepened density 
profile is sustained during the formation of protostar. Furthermore, the procedures for 
investigating the power-law index variation within a core suggest a steep index may present 
in the outside part of the core. 

One explanation for these results is that a core has to form a finite mass object before 
the inside-out collapse sets in. Vorobyov and Basu (2005, VB05) considered the effect of 
a finite mass reservoir during the core-forming and collapse phases. They found that an 
equivalent "rarefaction wave" is generated at the core boundary and propagates inward, 
which rapidly reduces the infall rate and steepens the density profile to p < — 2 (Figure 1 
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in VB05). Our results show a similar evolutionary trend with a steep power-law index seen 
in the outer region of our samples, as described in VB05. In addition, the VB05 model 
predicts that the infall velocity increases toward the center and then decreases when the 
velocity reaches a peak at ~ —0.2 to —0.4 km s _1 . A recent study on the velocity field 
of two starless cores show results consistent with the VB05 model (Figure 2 in Lee et al. 
2007). 

The physical reason for the "finite mass effect" could originate from the boundary of 
supercritical and subcritical region. For example, Tassis and Mouschovias (2007, TM07) 
show that during core formation via ambipolar diffusion, the magnetically supported mass 
cannot replenish the inner supercritical infall mass, resulting in a steep power-law index 
in the outer parts of the core (Figure 2d in TM07). Therefore, we may still see the cores 
smoothly blend in with the background but not confined by sharp edges as the literal 
meaning of finite mass. The origin of the "finite mass effect" and how significant is it for 
the star formation needs further investigation via observations on the density structure and 
also velocity structure of the cores at various evolutionary stages. 

We thank Shantanu Basu and Fred C. Adams for their insightful discussion on the 
first draft of this paper. This work is based on observations obtained at the Canada- 
France-Hawaii Telescope (CFHT) which is operated by the National Research Council 
of Canada, the Institut National des Sciences de l'Univers of the Centre National de la 
Recherche Scientifique of France, and the University of Hawaii. Access to the CFHT was 
made possible by the Institute of Astronomy and Astrophysics, Academia Sinica, Taiwan. 
CLH and SPL are supported by National Science Council of Taiwan under grant NSC 
96-2112-M-007-019-MY2 and NSC 98-21 12-M-007-007-MY3. 
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Fig. 1. — The extinction maps of the nine cores. The pseudocolor images display our CFHT 
H and Ks band and Spitzer's IRAC1 band data in blue, green, and red respectively. The 
contours represent the extinction Av and the contour levels are shown in the corner of each 
image. The black lines represent a the length scale of 0.1 pc. 
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Fig. 2. — The density profiles of the nine cores. The figures are the density profile of each 
core overlaid with the fitting results for the power-law density profile model in black solid 
lines. The numbers indicate the best-fit power-law index and its associated uncertainty. The 
dashed lines represent the uncertainties of best-fit profiles. The black dots show the density 
profile with the background extinction subtracted. The error bar is the standard deviation 
of annularly averaged extinction. 
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Fig. 3. — The best fit power-law index of the power-law density model. The black squares 
are the best-fit results, the black circles and the black triangles are obtained with fitting 
ranges ~ 0.2R — R and ~ 0.3R — R, respectively. 



